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We show how the AGBS model, originaUy developed for deep inelastic scattering applied to 
HERA data on the proton structure function, can also describe the RHIC data on single inclusive 
hadron yield for d + Au and p + p collisions through a new simultaneous fit. The single inclusive 
hadron production is modeled through the color glass condensate, which uses the quark(and gluon)- 
condensate amplitudes in momentum space. The AGBS model is also a momentum space model 
based on the asymptotic solutions of the BK equation, although a different definition of the Fourier 
transform is used. This aspect is overcome and a description entirely in transverse momentum of 
both processes arises for the first time. The small difference between the simultaneous fit and the 
one for HERA data alone suggests that the AGBS model describes very well both kind of processes 
and thus emerges as a good tool to investigate the inclusive hadron production data. We use this 
model for predictions at LHC energies, which agree very well with available experimental data. 
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I. INTRODUCTION 



The understanding of the high energy behavior of the quantum chromodynamics (QCD) scattering amplitudes 
towards saturation has increased over the last decades. The pioneering work of Gribov, Levin, and Ryskin 
which aimed keep the high energy scattering amplitudes unitary in the small- a; limit of QCD by means of parton 
recombination, has started the field of investigation of the nonlinear evolution equations of the QCD amplitudes. The 
called GLR equation has also settled the presence into its solutions of a semi-hard energy dependent scale Qs{x) 
that determines the onset of the saturation effects in the evolved amplitudes, which posteriorly was called saturation 
scale. For probes with transversal momentum kt < Qs, the saturation effects are important to unitarize the scattering 
amplitudes. 

Since the GLR equation, other nonlinear evolution equations were developed to describe this high energy limit of 
QCD [11,0], of which the simplest is the Balitsky-Kovchegov (BK) equation for the dipole-target amplitude. 

BK equation is seen as a "mean-field" approximation once it neglects higher order correlations between the parton 
content of the amplitudes that are present in the so-called Balitsky j^] and JIMWLK Q equations, which consider 
the target as a Color Glass Condensate (CGC). The BK equation describes the rapidity evolution of a color dipole 
scattering by a target, being the dipole formed by a quark-antiquark pair with the transverse coordinates x and y. 
In the fixed coupling case the scattering amplitude evolution in the two-dimensional position space reads 



where \xy\'^ = {x — y)^ is the dipole size, Y = Inl/x the rapidity variable and a = UsNc/tt. 

On the phenomenological side, the geometric scaling property Q observed in deep inelastic scattering (DIS) data 
have determined the presence on the cross sections of the rapidity dependent saturation scale Qs , as predicted from 
the GLR equation. This property establishes that the DIS cross sections depend on the Bjorken variable x and on the 
virtuality only through the combination x/Qs{x), where Qs is the saturation scale. An interesting feature between 
theory and phenomenology is the fact that only a saturation based model could describe the geometric scaling form 
of the DIS data at HERA [8] , giving us a strong empiric evidence of the saturation phenomena presence in the high 
energy collisions. 

Even though numerical analyses of the fixed and running coupling versions of the BK equation have been 
performed, giving us insight about the solutions, the analytical solutions could not be found. This situation has 
changed with the discovery of the equivalence fill] between the BK equation and the Fisher-Kolmogorov-Petrovsky- 
Piscounov (FKPP) equation |12j|. The later was largely studied in statistical physics [l^ [l3| and is known to admit 
traveling waves as asymptotic solutions, which translates into the geometric scaling property of the QCD when applied 
to BK equation. The analysis is not dependent either on the form of the nonlinear term of the equation neither on the 
form of the linear kernel, so that the property of the traveling wave solutions is said to be universal. Furthermore, at 
asymptotic rapidities, when the wave reaches its traveling behavior, it loses the memory about the form of the initial 
conditions. 

This procedure, called the traveling wave method of QCD, has received much theoretical effort, including studies 
of the higher order corrections of the BK kernel [l5l - [l7| . However, there are few applications to phenomenology. The 
goal of this work is to extend a traveling wave based model [31 to perform aiialyses both of the HERA DIS data on 
the proton structure function [l^ and of the RHIC data on hadron yield [13, through a simultaneous fit. In this 
way we can observe how the asymptotic solutions of the BK equation — as entering in the AGBS amplitude — behaves 
in nuclear collisions, with the parameters also limited to the DIS data at HERA. The result of this analysis is used 
to describe the LHC pt distribution data for p -f p collisions [l^, . 

The paper is organized as follows: in Section [H] we review the traveling wave solutions of the BK equation and its 
modeling to DIS phenomenology. Section IIIII is devoted to the mathematical formulation and kinematics applied to 
DIS at HERA and to hadron production at RHIC. Section ITVl shows the data set and the results of the simultaneous 
fit as well as its predictions for LHC and Section |V] presents our conclusions. 



II. THE TRAVELING WAVE SOLUTIONS OF BK EQUATION AND THE AGBS MODEL 

For the BK equation in the one-dimensional form, when the impact parameter is neglected, the amplitude in ([T]) 
takes the form Afyix, y) = Afyir), with r = \xy\. Thus it could be Fourier transformed through 



dYJ^Y{x,y) = a cPz 




- [Mvix^z) +AfYiz,y) ~ MY{x,y) ~ My{x, zjMviz.y)] , 



(1) 




— Mkr)U{r,Y), 



(2) 
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and the BK in momentum space reads [5| 

dyNik, Y) = axi~dL)Nik, Y) - aN{k, Y f, (3) 

where 

x(7)-27/;(l)-7/;(7)-V(l-7) (4) 

is the BFKL [g^l kernel and L = log(fc^/A;§), with a fixed soft scale fco. 

In the diffusive approximation for this kernel, i.e., taking the Taylor expansion of the kernel up to second order, 
it can be shown [llj that the BK equation ([3]) belongs to the same universality class as the FKPP equation, which 
means that the BK equation admits traveling waves N{L — VgaY) as solutions. The unique condition for the existence 
of these solutions is that the initial conditions N{L,Yo) decrease faster than exp(— 70^) for large L, with 70 > 7c- 
The color transparency property of QCD means that the scattering is weak when the dipole size is small, so that 
Af (X r"^. Using ([2]) to Fourier transform this equation it gives N cx k~'^'">, whit 70 = 1 at asymptotic L. The critical 
anomalous dimension 7c defines the group velocity Vg{'yc) of the asymptotic traveling wave. For the leading order 
BFKL kernel, 7c = 0.6275 so that the condition 70 > 7c is fulfilled. 

Through the traveling wave method it is obtained the asymptotic part of the BK amplitude, valid for large k. For 



the fixed coupling case it reads [11 



N{k,Y) 



log 



cxp 



W {kVQliY)) 



2ax"{lc)Y 



(5) 



where 



Of(F) = Q^exp( AF-— logF 



(6) 



is the saturation scale and A — avg — ax{lc)/^c- 

To make a phenomenological use of this asymptotic behavior, we need to add the infrared behavior for the amplitude 
in the region fc <C Qs. This was done in 18] using the Fourier transform of a two-dimensional step function, in order 
to unitarize the scattering amplitude in this saturated limit. The result is 



N{k,Y) (k) '"'^'^ c-log 



Qs{Y) 



(7) 



The equations ([S]) and ([7]) were analytically interpolated to form a momentum space parameterization of the dipole- 
target scattering amplitude, which reads 



N{k,Y) = 



1 / ^ , Qs 



(l-e-^-') 



(8) 



where 



Tdii = exp 



-7c log 



^rcd 



log' (2) 



2ax"hc)Y 



(9) 



with 



= log 1 



and 



QiiY) = k',e 



2 „XY 



(10) 



The equations above (I5 HTU1) describe the AGBS model [18] for the di pole -proton scattering amplitude, originally fitted 
to the HERA HI and ZEUS collaborations noncombined data [25l l26l|. The model was also used to investigate the 
effects of the fluctuation in the gluon number during the dipole evolution [27l - l30j at the HERA energies [3l| . 
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III. PHENOMENOLOGY WITH AGBS MODEL 
A. DIS at HERA 



A convenient frame to study the small- a; electron-proton deep inelastic scattering is the dipole frame, where the 
virtual photon emitted by the electron can split into a quark-antiquark pair that probes the target. In this way the 
photon-proton cross section factorizes in the photon wave function describing the probab ility of photon splitting into 
the qq pair and the dipole cross section describing the dipole interaction with the target [12, : 

all{Q\Y)= I d\ j\z\'fT,Lir,z;Q^)\^'Jdip{r,Y), (11) 

where z is the momentum fraction carried by the quark. Neglecting the impact parameter of this interaction, i.e., 
taking the target as a homogeneous disk with fixed radius, the dipole cross section turns out to be a function of the 
dipole-target forward amplitude. In the case of the e + p collisions we have 

adip(fc,r) =2^i?2^(fc,y), (12) 

where Rp is the radius of the target proton. 

The AGBS model was originally fitted to the proton structure function cx crj £ at HERA, which was written 
from (jlip in the momentum space as [isj 



)\'N{k,Y), (13) 



where N{k,Y) refers to the AGBS dipole amplitude and ^t,l = + to the photon wave function, written in 
momentum space through 

/^2 poo 7 

-^e'^ VvE'T,L(r,z)= / — Jo(fcr)vI/^,i(r, z). (14) 
(27r)^ Jq 27r 



Thus, we see from the convolution of the transforms of the photon wave function (|14p and of the dipole amplitude 
(p|) that the whole transform of the proton structure function has the form of a two-dimensional Fourier (or Hankel) 
transform. As we shall see, this fact will be important in our goal to perform a simultaneous fit of the AGBS amplitude 
to the HERA and RHIC data. 

B. Inclusive hadron production at RHIC 

Within the frame of the Color Glass Condensate (CGC), it was proven that the differential cross section for single- 
inclusive forward hadron production in high energy collisions depends only upon dipoles through j34l | 

dN K dz [ 2\ <-T fPt \ / 2\ 



dyh d'^pt (27r)2 J^^ z 



(15) 

+ Xlfg/p{Xl,Pt)NA (^'^2) Dh/g {z,Pt) 



where a summation over quark flavors q is understood. In this equation, pt and yh are the transverse momentum 
and rapidity of the produced hadron while fi/p and Dj^/^ refer to the parton distribution function of the incoming 
nucleon and to the hadron fragmentation function respectively, which are considered at the scale = p^ > 1 GeV^. 
Here we will use the CTEQ6 LO p.d.f's [sij and the LO KKP fragmentation functions [s^. For light hadrons the 
finite mass effects can be neglected, so that the pseudorapidity 77 and rapidity yh of the produced hadrons are similar 
77 « yh, giving the following kinematics: xp = y^ml+p'^ exp{r]h) / V Snn ~ Pt exp{yh)/VSNN, X2 = xiexp{-2yh) 
and xi = xp/z. 

The amplitude in momentum space is obtained through the Fourier transform of the dipole-proton forward scat- 
tering amplitude in the coordinate space 



iVf(^)(fc,y) = / dVe*AAj.(^)(r,r) =27r / dr r Mkr)^p^A^ir,Y), (16) 
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where Np and Na are amplitudes in the fundamental and adjoint representations describing, respectively, the scat- 
tering of the quark and giuon content inside the projectile wavefunction by the target. The Np amplitude for quarks 
is obtained from the adjoint Na by the replacement Ql — >■ {Cp/Ca)QI, with Cp/Ca — 4/9. 

In order to use the AGBS model in Eq. (fT5|) and fit it to RHIC data, we have to take into account the proper 
Fourier space where this equation was derived, given by Eq. ([T5)) . The AGBS amplitude N{k), calculated with a 
different definition of Fourier transform given by Eq. ([2]), can be translated to: 



iV(fc) = 27r 



(fN{k) 



1 dN{k) 
k dk 



(17) 



So, performing the first derivatives of the AGBS amplitude ([5]) we get analytically an amplitude in the appropriate 
Fourier space to describe the RHIC data on single inclusive hadron production through ([T5|) . 

Some conceptual comments about the fundamental and adjoint representations of the scattering amplitudes are in 
order. As pointed out in [s^, the dipole amplitude that enters in the calculation of DIS observables is an amplitude 
for quarks and so in the fundamental representation. Thus, in the dipole models applied to DIS, like the AGBS, the 
amplitude is a Np one without the Cp/Ca factor mentioned before. To avoid this disagreement between the DIS and 
heavy ion amplitudes in our simultaneous fit to HERA and RHIC data we use Eq. pT]) as a model for Np and scale 
Ql Ca/CpQI to obtain Na- 



IV. RESULTS 



HERA and RHIC data simultaneous fit 



Both HERA proton structure function F2 (with the amplitude (|S])) and RHIC hadron yield (|15p (using the 

amplitude ([TT])) data were used to simultaneously fit the AGBS model to DIS and heavy ion processes. 

In this analysis, all the last combined HERA data measurements of the proton structure function from HI and 
ZEUS Collaborations [l^ are fitted, within the following kinematic range: 



X < 0.01, 



(18) 



0.1 <Q^ < 150GcV^ 



(19) 



which corresponds to 244 data points. Both ranges include values of x low enough for the analysis to be in the high 
energy regime, and values of which allow us not to include DGLAP corrections, although higher values of virtuality 
do not change significantly the results. 

The RHIC data on single inclusive hadron production from BRAHMS 2l| and STAR [20] collaborations were also 
fitted. They were considered in the pt range higher than 1 GeV to guaranty that perturbative theory can be applied. 
To avoid contribution of large- a; in the target, we have focused our study in the forward (i/h > 2) rapidity region, 
although we also made the fit to the and mid-rapidity {yh ~ 1) region. This gives us a total of 22 data points in 
the forward region analysis and 38 data points when the mid-rapidity region was included. In order to get the equal 
contributions from the HERA data and from the RHIC data, we have assigned a weight of 11 (6 in the case including 
the region y^. — 1) to the latter. 

After all, for the simultaneous fit we have used 266 (only RHIC forward region) and 283 (RHIC forward and 
mid-rapidity regions) data points. Concerning the parameters, we keep fixed a — 0.2 and 7c — 0.6275, whose value 
corresponds to the one obtained from the LO BFKL kernel. The other parameters in the amplitude — A, fcg and 
x"i"fc) — are left to be free, as well as the proton radius Rp, which fixes the normalization of the dipole-proton cross 
section with respect to the dipole-proton amplitude for DIS at HERA, and the rapidity-dependent K factors. We 
have to fix one of the K factor to set the overall normalization; therefore, we set the value K{y}i = 4) = 0.7, which 
was obtained from two different dipole models applied to the RHIC data [H, . Only light quarks are considered 
and the values used for their masses were niu d s = 140 MeV. 



x7d-o.f. 


kl {xlO-3) 


A 


X"(7c) 


R(GcV-i) 


0.903 


1.13 ± 0.024 


0.165 ±0.002 


7.488 ± 0.081 


5.490 ± 0.039 



TABLE I. Parameters extracted from the fit to HI and ZEUS combined data [T^ on the proton structure function F2 at HERA. 



6 




123456 123456 
pJGeV] pJGeV] 



{a.)d + Au collisions at RHIC for yh > 1.0 {h)d + Au collisions at RHIC for yn > 2.2 



FIG. 1. Results for the RHIC charged hadron and tv yield for d + Au collisions from the simultaneous fit of AGBS to RHIC 
[13, HI and HERA ^ data. 





xVd.o.l. 


feo (xlO-3) 


A 


x"(7c) 


fl(GoV-i) 


K{y,, = 1.0) 


KivH = 2.2) 


KivH = 3.2) 


Kivh = 4.0) 


Vh > 2.2 


0.799 


2.760 ± 0.130 


0.190 ± 0.003 


5.285 ± 0.123 


4.174 ± 0.053 




2.816 ± 0.110 


2.390 ± 0.098 


0.7 


Vh > 1.0 


1.056 


1.660 ± 0.137 


0.186 ± 0.003 


6.698 ± 0.223 


4.695 ± 0.112 


6.172 ± 0.379 


3.783 ± 0.259 


3.256 ± 0.226 


0.7 



TABLE II. Parameters extracted from the simultaneous fit to HERA F2 (HI and ZEUS combined data [T^) and to the RHIC 
hadron yield for the d + Au collisions (BRAHMS and STAR data [13, Ell). The first line corresponds to the fit with the forward 
rapidity region of the RHIC data included, while in the second line the mid-rapidity region of the RHIC data set was also 
considered. 



First of all we performed a fit to DIS data, in order to verify the behavior of the AGBS model to the new combined 
HI and ZEUS data 19|. The values of the parameters, shown in the table[Tl are close to those obtained in the original 
AGBS fit to the HERA data, although the recent data set is a little bit different and does not need the inclusion of a 
normalization uncertainty of 5% in the HI data. This fit will also serve as a guideline for our simultaneous fit, once 
DIS processes are free from model dependent PDF's and fragmentation functions that can interfere on the best values 
for the model parameters. 

The simultaneous fit to DIS and hadron production in d + Au collisions has shown a good agreement with the 
data set, mainly when only forward RHIC data was considered, as seen in the values of the tabic HIl In the Figs. 



1(a) and 1(b) we see the description of the RHIC hadron yield for both lines of table [Hi The poorer description 



of the mid-rapidity region occurs because the target have not reached its gluon condensate state, when neither our 
amplitude nor the CGC formulation entering the Eq. are valid. The parameters values — even those got when the 
mid-rapidity data was included — shows that the AGBS model describes equally well the HERA and RHIC data, and 
the last improves the AGBS model. The saturation exponent A is slightly bigger than that coming from the HERA 
fit only, close to the expected value A = 0.2 0.3 demanded from NLO DIS [40|. The saturation scale is just a little 
bit greater than that extracted from the DIS fit only, which have given us an unexpected small value for the nuclear 
saturation scale. However, includes an additional factor A^f^ to account for the dipole interaction with the whole 
nuclear target and must be viewed as an average value. The value in the center of the nucleus should be higher and 
an analysis of the impact parameter influence on the AGBS amplitude will be important for the LHC energies, as 
seen from other saturation models [4l[. We have used A^s = 18.5 for d + Au collisions. 

Concerning the K factors, we had to fix the pion one to normalize such factors as explained before. We used 
the value K — 0.7 given by other two LO models [ssl. [39|. In fact, the first line of table HTl shows almost the same 
values of these LO models, with a factor of two for the hadron production {yh — 2.2 and yh = 3.2) as expected from 
LO models. We must note, however, that these factors are just normalization factors over the equation (|15p in the 
fitting procedure, describing the model uncertainties both on the DIS and hadron yield data, that actualy has 15% 
of normalization uncertainty. 

The results for the simultaneous fit considering p + p collisions at RHIC are presented in the Figs. 2(a) and |2(b)] 
where we see an even worse description in the fit of the mid-rapidity data in comparison with the d+Au collisions. This 
is because the gold nucleus target has a greater partonic content than the proton at the same rapidity, being better 
described with the CGC formulation used here. The same does not happen with the proton target and corrections 
to Eq. (fT5|) should be important to describe this data at mid-rapidities. All in all, considering this fact and also 
remembering that AGBS is a model for the small-x region of the QCD, the fit to the forward rapidities is the best 
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choice for future use of the AGBS model in the single inclusive hadron production. 



BRAHMS data x200, (h^ + h")/2, yu=1 .0 ^ 

BRAHMS data xlOO, h", yh=2.2 ^ 

BRAHMS data x50, K, ^^,=3.2 i 

STAR data x20, yh=3.3 t 

STAR data x1 D, Ti„, yh=3.8 t 

STAR data, ii , yh=4.0 ^ 




p,[GeV] 




BRAHMS dataxlOO, ti", 



STAR data x10,<, 



STAR data. Jt , yh=4.0 ^ 



p,[GeV] 



(a)p + p collisions at RHIC for j/h > 1.0 



(b)p + p collisions at RHIC for > 2.2 



FIG. 2. Results for the RHIC charged hadron and tt" yield for p+p collisions from the simultaneous fit of AGBS to RHIC and 
HERA data. 



B. LHC predictions 



In this section we present predictions of the AGBS model ioi p + p and p + Pb collisions at LHC energies of 14 and 
8.8 TeV, respectively. We use the parameters obtained from the simultaneous fit in the case of forward rapidity region 
of RHIC data (first line of the table HII as explained above and our results are shown in the Figs. 3(a) and |3(b)] 

Fig. m shows the description of the recent LHC CMS [l^, [1^ pt distribution data for proton-proton collisions using 
our model. We used the following expression to change the variable from rapidity to the measured pseudorapidity: 



y{v,Pt,m) = 2 



ni? + cosh 7y + pt sinh i] 



m? + pI cosh rj — pt sinh rj 



(20) 



where m — 0.139 GcV is the pion mass. We observed that for light hadrons (m up to 0.4 GeV), as in the case of the 
RHIC data, using rj k, does not affect considerably the results. 

Surprisingly, a very good description was obtained, although the K factors are large. This could be explained as 
an uncertainty of the AGBS model in the comparison with the pseudorapidity averaged data performed by the CMS 
collaboration. As the AGBS model is not supposed to describe the central rapidity region since it is a model for the 
low- a: part of the target (projectile), we used a averaged value 77 = 1.4 so that the description of the pseudorapidity 
averaged data over the region \r]\ < 2.4 should imply some disagreements, which are in the K factor. 

Standing for NLO corrections, the K factors account for the qq and gq interactions, while Eq. ((T5]) considers only 
gluons in the target. Thus, as expected, the K values decrease with the energy, because the target (projectile) gluon 
wavefunction is more important when compared with the quark density as the energy increases. In other words, Eq. 
(1151) should receive less quark corrections in the target for higher energies. 



V. CONCLUSIONS AND DISCUSSION 



We have shown in this work the compatibility of the AGBS model with both the single inclusive hadron yield 
at RHIC and the small-x DIS at HERA. Through Eq. ((T7)) one could write analytically the AGBS model in the 
appropriate Fourier space used to describe the single inclusive hadron yield from the Color Glass Condensate. We 
have performed a new simultaneous fit of the AGBS model to the F2 at HERA and hadron yield at RHIC, which 
results agreed with the one performed to HERA data alone, meaning that our expression correctly describes the AGBS 
in the standard Fourier space of the CGC formalism applied to hadron production. 

The model describes well the forward rapidity region of the RHIC data for d + Au and p + p collisions, but fails 
at mid-rapidities. Thus, being a model for the small- a; dipole amplitude, the parameters from the fit for the forward 
rapidities should be used to further employ of the AGBS model in inclusive hadron production. The same formulation 
was used to make predictions to the LHC p + p and p + Pb collisions at 14 and 8.8 TeV, respectively. As a matter of 
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pJGeV] 

(a)p + p collisions at LHC 




pJGaV] 

(b)p + Pb collisions at LHC 



FIG. 3. Predictions of the AGBS model to the LHC charged hadron yield for p + p and p + Pb collisions. Parameters from the 
simultaneous fit of AGBS to forward rapidity RHIC and HERA data (first line of table |TT|. We have used AaB = 20 for lead 
target. 



7TeVx4, K=8.5 h 
2.36 TeV x 2, K=9.7 
0.9 TeV, K=11 




FIG. 4. Predictions of the AGBS model to the LHC CMS charged hadron yield for p + p collisions at a/s = 0.9, 2.36 and 7 
TeV. The experimental points are from CMS for I77I < 2.4 [2^. [23l]. Parameters extracted from the simultaneous fit of AGBS 
to forward rapidity RHIC and HERA data (first line of table HI)). 



data comparison, the model was also used to describe the CMS data on the hadron yield ior p + p collisions [22, |23| 
at different energies. Fig. U shows a pretty good description of the data, despite the high values of the K factors. 
As explained, this could be viewed as uncertainties of the model for the forward amplitude in the central region of 
pseudorapidity since the CMS data sums the entire region I77I < 2.4, and also in the target quark corrections not 
included in the CGC formulation for the hadron production. 

Although a good description of the nuclear data on inclusive hadron production was obtained, a study of the impact 
parameter dependence of the AGBS dipolc amplitude is worth to investigate how sensitive in the transverse plane is 
the model. A better understanding of the parameters now averaged over the impact parameter would improve the 
model and its applicability to the high multiplicity events at LHC [4l| . 

Furthermore, in order to get rid of the hadron fragmentation effects in our calculations, it is interesting to investigate 
prompt photon production in hadron collisions. As they interact only electromagnetically with the medium, the direct 
(prompt) photon arises as a good tool to investigate initial state effects and saturation effects in hadronic collisions. 
This observable is well studied in the color dipole formalism [42] and could be also used to investigate the heavy quark 
contribution to the dipole scatterin g am plitudes through the photon plus heavy quark production, as recently done 
in [il] with both dipole models [ssl. l39j. 

At last, a next-to-leading order (NLO) phenomenological modeling of the traveling wave QCD method is needed 
to mark the region of applicability of the AGBS model. The theory of the method at higher orders has already been 
done, in the spirit of regularization group scheme [Tsl - fTTj . However there are no applications in phenomenology yet, 
being an interesting field for future studies. 
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